Dataset

The data for following exercises stem from the publication Grassland ecosystem recovery after soil disturbance depends on nutrient supply rate and are publicly available at Dryad. The data were obtained during the long-term field experiment Cedar Creek LTER and target the effects of human disturbances on grassland ecosystem functioning and biodiversity.

Data collection

Treatment

The disturbance treatment was replicated in the three old-fields (A, B and C) in a completely randomised block design (two treatments in each of three fields for a total of 6 35 × 55 m large plots). In 1982, in each of the fields, one of these two 35 × 55 m areas was selected to be disturbed with a disk harrow. After, the soil was hand raked to smooth the soil and remove any remaining vegetation, so that subsequent colonisation was solely from seeds or small rhizome fragments. Within each of the 6 large plots, the 54 small plots were arrayed in 6 × 9 grid with 1 m buffers between each plot. Aluminium flashing was buried to depth of 30 cm around each plot to prevent horizontal movement of nutrients and spreading of plants through vegetative growth.

The nutrient treatments were replicated six times in a completely randomised design in each of the 35 × 55 m plots (54 4 × 4 m small plots) yielding 324 (6 x 54) plots. The analyses focuses on two nutrient treatments:

  1. Control (no nutrients; Treatment I) and
  2. Other Nutrients and 9.5 g of N (Treatment F)

Sampling and analysis

At peak biomass (mid-July to late August), all aboveground biomass was clipped in a 3 m by 10 cm strip (0.3 m2) in each plot. Note that there were 4 years when the disturbed plots were not sampled or only sampled in a single field. The biomass was sorted into dead, previous year’s growth (litter) and current year’s growth (live biomass). Live biomass was sorted to species, dried and weighed. We estimated total aboveground biomass as the summed biomass of all non-woody species in each 0.3 m2 sample, converted to g/m2.

Species richness is the number of species in each 0.3 m2 sample. We quantified plant diversity as the Effective Number of Species based on the Probability of Interspecific Encounter (ENSPIE), a measure of diversity that is more robust to the effects of sampling scale and less sensitive to the presence of rare species than species richness. ENSPIE is equivalent to the Inverse Simpson’s index of diversity (\(1 / \sum_{i=1}^{S} p_i^2\) where \(S\) is the total number of species and \(p_i\) is the proportion of the community biomass reesented by species \(i\)).

Load data

seabloom <- read.table(here("2_Modeling/Data_preparation/seabloom-2020-ele-dryad-data/cdr-e001-e002-output-data.csv"),
                       sep = ",", header = TRUE)

Explore data

dim(seabloom)
## [1] 5040   16
str(seabloom)
## 'data.frame':    5040 obs. of  16 variables:
##  $ exp       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ field     : chr  "A" "A" "A" "A" ...
##  $ plot      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ disk      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ yr.plowed : int  1968 1968 1968 1968 1968 1968 1968 1968 1968 1968 ...
##  $ ntrt      : int  9 9 9 9 9 9 9 9 9 9 ...
##  $ nadd      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ other.add : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ year      : int  1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 ...
##  $ dur       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ precip.mm : num  879 858 876 994 859 ...
##  $ precip.gs : num  374 551 527 503 512 ...
##  $ mass.above: num  61.3 135.9 216.9 302.8 586.9 ...
##  $ rich      : int  13 15 14 14 16 14 8 7 9 13 ...
##  $ even      : num  0.463 0.156 0.204 0.14 0.269 ...
##  $ ens.pie   : num  6.02 2.34 2.85 1.96 4.31 ...

Exploring the data reveals 16 variables with each 5040 data points:

  • exp: treatments in split-plot design: 1 = disturbance (Control or Disked, 35 × 55 m plots) and 2 = nutrient addition (9 levels, 4 × 4 m plots)
  • field: three experimental fields A, B and C
  • plot: 54 plots within fields
  • disk: disking treatment (0 = intact at start of experiment, 1 = disked at start of experiment)
  • yr.plowed: last year field was plowed for agriculture (A: 1968, B: 1957 and C: 1934)
  • ntrt: nine levels representing different combinations of nitrogen (0 to 27.2 g N year-1 added as NH4NO3) and other nutrients (20 g m−2 year−1 P205; 20 g m−2 year−1 K20; 40 g m−2 year−1 CaCO3; 30.0 g m−2 year−1 MgSO4; 18 μg m−2 year−1 CuSO4; 37.7 μg m−2 year−1 ZnSO4; 15.3 μg m−2 year−1 CoCO2; 322 μg m−2 year−1 MnCl2 and 15.1 μg m−2 year−1 NaMoO4; details see Table S1 in publication). Nutrients were applied twice per year in mid-May and mid-June.
  • nadd: nitrogen additon rate (g/m2/yr)
  • other.add: other nutrient treatment (0 = control, 1 = other nutrients added)
  • year: sampling year
  • dur: duration of experiment
  • precip.mm: annual precipitation (mm)
  • precip.gs: growing season precipitation (mm)
  • mass.above: aboveground biomass (g/m2)
  • rich: species richness (species/0.3 m2)
  • even: Simpson’s evenness
  • ens.pie: effective number of species, (= probability of interspecific encounter decimal equivalent to inverse Simpson’s diversity)

Overview

The pairs function yields an overview over the numerical data that we will use for the following exercises.

pairs(seabloom[, -c(1:3, 5:6, 8:10, 12, 16)],
      lower.panel = NULL)

Metamodel

A metamodel summarizes the concept behind a model and links it to theory. Here, the metamodel is visualized as a directed acyclic graph (DAG) which reads as: productivity (biomass) is directly influenced by the environment (nutrients, disturbance and precipitation) on the one hand and biodiversity (richness and evenness) on the other hand. Also some elements of the environment influence biodiversity and thus, have an indirect effect on productivity via biodiversity.

Subset to one year

For simplicity, set the focus on only one year. Then, 240 observations remain.

seabloom <- seabloom[seabloom$year == 2000, ]
dim(seabloom)
## [1] 240  16

Linear model

First, implement the metamodel into a linear model (LM). For this, three models are necessary: one that accounts for the direct- and two for the indirect effects.

Direct effects

If a linear model is visualized in a DAG, it becomes apparent that there is a set of permitted, but unanalyzed correlations among the predictors. These correlations have a huge influence on the coefficients between the predictors \(x_{1\ldots i}\) (here, nutrients, precipitation, richness and evenness) and the response \(y\) (here, biomass). Despite their importance, it is impossible to include any information about the reason of the correlations between the predictors. Further, they make it nearly impossible to create a proper causal model, since there are many possible causal relations that can create a set of “unanalyzed associations” what hinders interpretations (James B. Grace 2021).

Exercise: write the linear model in the DAG above in R syntax.

# lm.dir <- lm(mass.above ~ nadd + precip.mm + rich + even,
#              data = seabloom)
# summary(lm.dir)

We encounter a problem: the estimates for precip.mm are all NAs. Reason is the focus on a single year (i.e. the year 2000) what fixes the value for precipitation.mm to 599.7 mm (i.e. the precipitation in 2000) and as it is impossible to predict anything from a single value, NA is returned.

summary(seabloom$precip.mm)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   599.7   599.7   599.7   599.7   599.7   599.7

Thus, we exclude precipitation.mm from the model from now on.

lm.dir <- lm(mass.above ~ nadd + rich + even, data = seabloom)
summary(lm.dir)
## 
## Call:
## lm(formula = mass.above ~ nadd + rich + even, data = seabloom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -218.85  -59.09  -15.35   36.31  689.30 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  73.4044    38.5707   1.903 0.058242 .  
## nadd          6.7044     0.8345   8.034  4.5e-14 ***
## rich         12.7405     3.6505   3.490 0.000576 ***
## even         36.8867    40.5034   0.911 0.363379    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 100.4 on 236 degrees of freedom
## Multiple R-squared:  0.2153, Adjusted R-squared:  0.2053 
## F-statistic: 21.58 on 3 and 236 DF,  p-value: 2.187e-12

Indirect effects

To account for the indirect effects, two additional LMs are necessary: one with richness and one with evenness as response.

lm.rich <- lm(rich ~ nadd, data = seabloom)
summary(lm.rich)
## 
## Call:
## lm(formula = rich ~ nadd, data = seabloom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0318 -1.7746 -0.5349  1.2254  8.9682 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.03177    0.19477  30.969  < 2e-16 ***
## nadd        -0.12856    0.01614  -7.965 6.81e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.216 on 238 degrees of freedom
## Multiple R-squared:  0.2105, Adjusted R-squared:  0.2071 
## F-statistic: 63.44 on 1 and 238 DF,  p-value: 6.81e-14
lm.even <- lm(even ~ nadd, data = seabloom)
summary(lm.even)
## 
## Call:
## lm(formula = even ~ nadd, data = seabloom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39021 -0.13808 -0.03255  0.11926  0.50176 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.494342   0.017554  28.161   <2e-16 ***
## nadd        0.003422   0.001455   2.353   0.0195 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1998 on 238 degrees of freedom
## Multiple R-squared:  0.02273,    Adjusted R-squared:  0.01862 
## F-statistic: 5.534 on 1 and 238 DF,  p-value: 0.01946

Conclusion

The direct effect model showed that in the year 2000 biomass was statistically significantly positively related to nutrient input (the nitrogen additon rate nadd) and species richness (rich), but not to evenness (even). Additionally, nutrient input had a statistically significantly negative influence on species richness and to a lesser extent, also on evenness.

Structural equation modeling

Structural equations aspire to represent cause-effect relationships and thus, they can be used to represent scientific, causal hypotheses. A key feature of structural equation models is the ability to investigate the networks of connections among system components (James B. Grace et al. 2012).

To evaluate the SEMs, we will use the package lavaan. This powerful piece of software relies on the computation of covariance matrices to fit the structural equations. It comes with full support for categorical data (any mixture of binary, ordinal and continuous observed variables), can handle both latent and composite variables, performs mediation analysis and calculates the magnitude of indirect effects (Rosseel 2012, 2021).

library("lavaan")
## This is lavaan 0.6-9
## lavaan is FREE software! Please report any bugs.

Collinearity

Variables that are highly correlated with \(r > 0.85\) become redundant. Then, one of them can be dropped or they can be modeled summarized as a latent variable (James B. Grace 2006).

Inspecting our set of variables shows that correlations between all variables are fairly low.

round(cor(seabloom[, c(4, 7, 13:15)]), digits = 2)
##             disk  nadd mass.above  rich  even
## disk        1.00  0.00       0.05 -0.01  0.02
## nadd        0.00  1.00       0.41 -0.46  0.15
## mass.above  0.05  0.41       1.00  0.00 -0.02
## rich       -0.01 -0.46       0.00  1.00 -0.59
## even        0.02  0.15      -0.02 -0.59  1.00

An example SEM

This toy model shall illustrate the logic of SEM. It contains the same variables and aims at evaluating the same metamodel as the LM before to ease Comparability. A huge benefit of SEM is that variables can appear as both predictors and responses what allows the evaluation of direct and indirect effects in one go (as shown in the directed acyclic graph (DAG) below). In SEM jargon, predictors (nodes that have no arrows pointing at them) are colled exogenous-, while responses (nodes that have arrows pointing at them) endogeneous variables.

The table below summarizes the operators of lavaan syntax. For example, an arrow in a DAG is represented by a tilde.

Formula type Operator Meaning Example
regression ~ is regressed on y ~ x
correlation ~~ correlate errors for y1 ~~ y2
latent variable =~ set reflective indicators Height =~ y1 + y2 + y3
composite variable <~ set formative indicators Comp1 <~ 1*x1 + x2 + x3
mean/intercept ~ 1 estimate mean for y without xs y ~ 1
label parameter * name coefficients y ~ b1*x1 + b2*x2
define quantity := define quantity TotalEffect := b1*b3 + b2
define thresholds | for ordered categorial variable u | t1 + t2

Exercise: translate the model in the DAG above into lavaan syntax with the help of this table.

simple <-
"mass.above ~ nadd + disk + rich + even
rich ~ nadd
even ~ nadd"

Normal distribution

lavaan uses a \(\chi^2\) test to compare the estimated- to the observed covariance matrix to compute the goodness of fit for the SE model under the assumption that all observations are independent and all variables follow a (multivariate) normal distribution (James B. Grace 2006). Note, that these distributional assumptions only apply to endogenous variables, whereas the distribution of exogenous variables has no bearing on assumptions (James B. Grace 2021).

We use a graphical method to assess the fit of the endogeneous variables to a normal distribution, the quantile-quantile plots (Q-Q plots). Hereby, the quantiles of the data are compared to those of a theoretical distribution (i.e., the normal distribution). If the data would be normally distributed, the points would match one to one and thus, align diagonally. In the Q-Q plot, this match is indicated by the black line.

The package MVN allows to plot several Q-Q plots at once and also offers several tests to multivariate normal distribution (MVN). Here, we employ the Henze-Zirkler’s test that has been recommended as a formal test of MVN (Mecklin and Mundfrom 2005).

library("MVN")

mvn(data = seabloom[, c(13:15)], mvnTest = "hz", univariatePlot = "qqplot")
## $multivariateNormality
##            Test       HZ p value MVN
## 1 Henze-Zirkler 4.527585       0  NO
## 
## $univariateNormality
##               Test   Variable Statistic   p value Normality
## 1 Anderson-Darling mass.above    6.5613  <0.001      NO    
## 2 Anderson-Darling    rich       4.8536  <0.001      NO    
## 3 Anderson-Darling    even       3.5399  <0.001      NO    
## 
## $Descriptives
##              n       Mean     Std.Dev      Median       Min     Max        25th
## mass.above 240 211.002563 112.5962210 179.7370000 26.966000 897.564 140.1585000
## rich       240   4.979167   2.4892201   5.0000000  1.000000  15.000   3.0000000
## even       240   0.522363   0.2016525   0.4962043  0.197217   1.000   0.3741316
##                  75th      Skew   Kurtosis
## mass.above 265.128750 2.1849256  8.7573782
## rich         6.000000 0.9118942  0.7982559
## even         0.636105 0.7113656 -0.1130884

Alternatively, also histograms in comparison overlaid with a normal distribution allows to gain insight into the distribution:

par(mfrow = c(1, 3))
hist(seabloom$mass.above, prob = TRUE, main = "")
x <- seq(min(seabloom$mass.above), max(seabloom$mass.above), length = 400)
f <- dnorm(x, mean = mean(seabloom$mass.above), sd = sd(seabloom$mass.above))
lines(x, f, col = "red", lwd = 2)


hist(seabloom$rich, prob = TRUE, main = "")
x <- seq(min(seabloom$rich), max(seabloom$rich), length = 400)
f <- dnorm(x, mean = mean(seabloom$rich), sd = sd(seabloom$rich))
lines(x, f, col = "red", lwd = 2)


hist(seabloom$even, prob = TRUE, main = "")
x <- seq(min(seabloom$even), max(seabloom$even), length = 400)
f <- dnorm(x, mean = mean(seabloom$even), sd = sd(seabloom$even))
lines(x, f, col = "red", lwd = 2)

Exercise: would you infer that the endogenous variables meet the assumption of being (multivariate) normally distributed from the results of the the Q-Q plots, the Henze-Zirkler test and the histograms? And why (not)?

Data transformation

Transformation of the three endogenous variables shall make them more “normal.” However, this works only for continuous variables. Count data (in our example rich) are discrete data and thus, only have non-zero probabilities on the natural numbers.

par(mfrow = c(1, 3))

# log transform for the continuous variables
seabloom$log.mass.above <- log(seabloom$mass.above)
seabloom$log.even <- log(seabloom$even)
seabloom$log.rich <- log(seabloom$rich)


par(mfrow = c(1, 3))
hist(seabloom$log.mass.above, prob = TRUE, main = "")
x <- seq(min(seabloom$log.mass.above), max(seabloom$log.mass.above), length = 400)
f <- dnorm(x, mean = mean(seabloom$log.mass.above),
           sd = sd(seabloom$log.mass.above))
lines(x, f, col = "red", lwd = 2)


hist(seabloom$log.rich, prob = TRUE, main = "")
x <- seq(min(seabloom$log.rich), max(seabloom$log.rich), length = 400)
f <- dnorm(x, mean = mean(seabloom$log.rich), sd = sd(seabloom$log.rich))
lines(x, f, col = "red", lwd = 2)


hist(seabloom$log.even, prob = TRUE, main = "")
x <- seq(min(seabloom$log.even), max(seabloom$log.even), length = 400)
f <- dnorm(x, mean = mean(seabloom$log.even), sd = sd(seabloom$log.even))
lines(x, f, col = "red", lwd = 2)

mvn(data = seabloom[, c(17:19)], mvnTest = "hz", univariatePlot = "qqplot")
## $multivariateNormality
##            Test       HZ      p value MVN
## 1 Henze-Zirkler 1.612818 1.866308e-05  NO
## 
## $univariateNormality
##               Test       Variable Statistic   p value Normality
## 1 Anderson-Darling log.mass.above    0.7093  0.0633      YES   
## 2 Anderson-Darling    log.even       0.5155  0.1893      YES   
## 3 Anderson-Darling    log.rich       3.8487  <0.001      NO    
## 
## $Descriptives
##                  n       Mean   Std.Dev     Median       Min      Max
## log.mass.above 240  5.2309769 0.4976908  5.1914945  3.294577 6.799684
## log.even       240 -0.7231284 0.3887750 -0.7007676 -1.623451 0.000000
## log.rich       240  1.4787594 0.5167328  1.6094379  0.000000 2.708050
##                      25th       75th       Skew   Kurtosis
## log.mass.above  4.9427706  5.5801888 -0.2658087  1.3687570
## log.even       -0.9831481 -0.4523943 -0.1125629 -0.4850788
## log.rich        1.0986123  1.7917595 -0.2279647 -0.5296575

Fit the model

Now, let’s fit the model with lavaan’s sem function.

fit.simple <- sem(simple, data = seabloom)
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan
## WARNING: some observed variances are (at least) a factor 1000 times larger than
## others; use varTable(fit) to investigate

Oups, the algorithm converges with a warning. Kindly, it informs us how to fix this.

Exercise: let’s obey the software and execute the code from the hint.

# varTable(fit.simple)

This reveals an enormous difference in magnitude between the variance of biomass (mass.above) and the other variables.

Rescale variables

To remove this difference in magnitude between the variables, we divide mass.above by 100 what changes the unit from g/m2 to 10 mg/m2. The boxplots show that the range of the variables is now more similar.

seabloom$mass.above <- seabloom$mass.above / 100
boxplot(seabloom[, c(4, 7, 13:15)])

Then, run sem again with the rescaled biomass variable:

fit.simple <- sem(simple, data = seabloom)
summary(fit.simple, fit.measures = TRUE)
## lavaan 0.6-9 ended normally after 39 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         9
##                                                       
##   Number of observations                           240
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               104.782
##   Degrees of freedom                                 3
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               226.074
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.531
##   Tucker-Lewis Index (TLI)                      -0.407
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -822.542
##   Loglikelihood unrestricted model (H1)       -770.151
##                                                       
##   Akaike (AIC)                                1663.085
##   Bayesian (BIC)                              1694.411
##   Sample-size adjusted Bayesian (BIC)         1665.883
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.376
##   90 Percent confidence interval - lower         0.316
##   90 Percent confidence interval - upper         0.439
##   P-value RMSEA <= 0.05                          0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.141
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   mass.above ~                                        
##     nadd              0.067    0.008    8.161    0.000
##     disk              0.123    0.131    0.941    0.346
##     rich              0.127    0.029    4.381    0.000
##     even              0.361    0.322    1.119    0.263
##   rich ~                                              
##     nadd             -0.129    0.016   -7.998    0.000
##   even ~                                              
##     nadd              0.003    0.001    2.362    0.018
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .mass.above        0.987    0.090   10.954    0.000
##    .rich              4.872    0.445   10.954    0.000
##    .even              0.040    0.004   10.954    0.000

Goodness of fit

This time, the model converged, however, with poor fit:

  • The ratio of the test statistic and the degrees of freedom should be smaller than 2. Here, the ratio is 104.78 / 3 = 34.93, what indicates that the model is quite far away from a decent fit (James B. Grace 2021).
  • The \(p\)-value, which represents the probability of the data given our model (or no significant discrepancy between model and data), should be larger than 0.05 (James B. Grace 2021).
  • The comparative fit index (CFI) ranges between zero and one, whereas a value \(> .95\) is considered as a good fit indicator (Hu and Bentler 1999). Our CFI of 0.54 is fairly below this threshold.

Modification indices

To improve the model fit, we look for missing paths via modification indices. They indicate an estimated drop in the model \(\chi^2\) resulting from freeing fixed parameters via including a missing path. 3.84 is considered as the critical threshold, the “single-degree-of-freedom \(\chi^2\) criterion” (James B. Grace 2021).

modindices(fit.simple, minimum.value = 3.84)
##     lhs op        rhs     mi     epc sepc.lv sepc.all sepc.nox
## 15 rich ~~       even 84.811  -0.261  -0.261   -0.594   -0.594
## 16 rich  ~ mass.above 51.301 -10.881 -10.881   -4.969   -4.969
## 17 rich  ~       even 84.811  -6.596  -6.596   -0.534   -0.534
## 19 even  ~ mass.above 79.659  -0.399  -0.399   -2.248   -2.248
## 20 even  ~       rich 84.811  -0.054  -0.054   -0.661   -0.661

In the column mi (for modification index) we look for high values. Note, however, that the modification indices are uninformed suggestions and further adaptations of the model based on their information needs to be based on theory.

In this example, the modification indices indicate–amongst other–a missing relation between richness and evenness. Including this relation into the model would improve its fit by a change in the \(\chi^2\) by 84.811.

This path is necessary as richness and evenness are computationally related to each other (they are not independent quantities). Thus, let’s include a correlation between rich and even (Note: in SEM, a correlation between two variables points to an omitted common cause/variable that drives this correlation). With the function update() it is possible to directly incorporate the missing path into the specified model without rewriting it from scratch:

fit.simple.up <- update(fit.simple, add = "rich ~~ even")
summary(fit.simple.up, fit.measures = TRUE)
## lavaan 0.6-9 ended normally after 45 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        10
##                                                       
##   Number of observations                           240
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 0.143
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.931
## 
## Model Test Baseline Model:
## 
##   Test statistic                               226.074
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000
##   Tucker-Lewis Index (TLI)                       1.039
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -770.223
##   Loglikelihood unrestricted model (H1)       -770.151
##                                                       
##   Akaike (AIC)                                1560.445
##   Bayesian (BIC)                              1595.252
##   Sample-size adjusted Bayesian (BIC)         1563.554
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.037
##   P-value RMSEA <= 0.05                          0.961
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.007
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   mass.above ~                                        
##     nadd              0.067    0.008    8.118    0.000
##     disk              0.123    0.131    0.941    0.346
##     rich              0.127    0.036    3.523    0.000
##     even              0.361    0.401    0.900    0.368
##   rich ~                                              
##     nadd             -0.129    0.016   -7.998    0.000
##   even ~                                              
##     nadd              0.003    0.001    2.362    0.018
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##  .rich ~~                                             
##    .even             -0.261    0.033   -7.916    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .mass.above        0.987    0.090   10.954    0.000
##    .rich              4.872    0.445   10.954    0.000
##    .even              0.040    0.004   10.954    0.000
modindices(fit.simple.up, minimum.value = 3)
## [1] lhs      op       rhs      mi       epc      sepc.lv  sepc.all sepc.nox
## <0 rows> (or 0-length row.names)

Now, the model has a decent fit with a ratio of test statistic and degrees of freedom much smaller than two (i.e. 0.07) and a \(p\)-value of 0.93. Further, also the CFI is now 1.

Another look at the modification indices shows that further modifications would yield only smallest improvements to the model fit, so that we can ignore them confidently:

modindices(fit.simple.up, minimum.value = 0)
##     lhs op        rhs    mi    epc sepc.lv sepc.all sepc.nox
## 12 nadd ~~       disk 0.000  0.000   0.000       NA    0.000
## 16 rich  ~ mass.above 0.003  0.102   0.102    0.046    0.046
## 18 rich  ~       disk 0.003  0.013   0.013    0.002    0.005
## 19 even  ~ mass.above 0.111  0.057   0.057    0.318    0.318
## 21 even  ~       disk 0.111  0.007   0.007    0.017    0.035
## 22 nadd  ~ mass.above 0.000  0.000   0.000    0.000    0.000
## 25 nadd  ~       disk 0.000  0.000   0.000    0.000    0.000
## 26 disk  ~ mass.above 0.585 -0.789  -0.789   -1.809   -1.809
## 27 disk  ~       rich 0.020 -0.002  -0.002   -0.008   -0.008
## 28 disk  ~       even 0.133  0.056   0.056    0.023    0.023
## 29 disk  ~       nadd 0.000  0.000   0.000    0.000    0.000

Standardized coefficients

Using analysis of covariances allows for estimation of both unstandardized (raw) and standardized coefficients. While the analysis is based on covariances for estimating unstandardized coefficients, it is based on correlations for estimating standardized coefficients. The computational relation between correlations and covariances is:

\(r_{xy} = \frac{cov_{xy}}{SD_x \times SD_y}\)

… where \(r_{xy}\) are the correlations, \(cov_{xy}\) the covariances and \(SD_x\) and \(SD_y\) the respective standard deviations (the square roots of the respective variance \(= \sqrt{var_x}\) and \(= \sqrt{var_y}\)) (James B. Grace 2021).

Significance testing of standardized coefficients violates statistical principles as significance tests are based on unstandardized (raw) coefficients Thus, standardized coefficients are used for interpretation only (James B. Grace 2021).

For raw coefficients, the predicted effects are in raw units and have a pretty straight-forward interpretation: as in our example, e.g., the predicted change in 10 mg/m2 above-ground biomass associated with a certain change in added nutrients. For the standardized coefficients, however, the interpretation is more complex: they express the predicted changes in terms of standard deviation units. Thus, they are only interpretable within a sample what hinders generalization. Still, they ease comparison as they are in the same units across various pathways (James B. Grace 2021).

There are different standardizations available in lavaan, whereof std.all is based on the standard deviation:

standardizedsolution(fit.simple.up, type = "std.all")
##           lhs op        rhs est.std    se       z pvalue ci.lower ci.upper
## 1  mass.above  ~       nadd   0.529 0.055   9.676  0.000    0.422    0.636
## 2  mass.above  ~       disk   0.054 0.057   0.943  0.346   -0.058    0.165
## 3  mass.above  ~       rich   0.281 0.079   3.575  0.000    0.127    0.436
## 4  mass.above  ~       even   0.065 0.072   0.901  0.368   -0.076    0.205
## 5        rich  ~       nadd  -0.459 0.048  -9.516  0.000   -0.553   -0.364
## 6        even  ~       nadd   0.151 0.063   2.403  0.016    0.028    0.274
## 7        rich ~~       even  -0.594 0.042 -14.242  0.000   -0.676   -0.513
## 8  mass.above ~~ mass.above   0.782 0.046  17.168  0.000    0.692    0.871
## 9        rich ~~       rich   0.790 0.044  17.850  0.000    0.703    0.876
## 10       even ~~       even   0.977 0.019  51.677  0.000    0.940    1.014
## 11       nadd ~~       nadd   1.000 0.000      NA     NA    1.000    1.000
## 12       nadd ~~       disk   0.000 0.000      NA     NA    0.000    0.000
## 13       disk ~~       disk   1.000 0.000      NA     NA    1.000    1.000

Model comparison

To evaluate whether the simple or the updated model perform better, we can calculate the Akaike information criterion (AIC) and compute a significance test based on the \(\chi^2\) test statistic with the anova function.

anova(fit.simple, fit.simple.up)
## Chi-Squared Difference Test
## 
##               Df    AIC    BIC    Chisq Chisq diff Df diff Pr(>Chisq)    
## fit.simple.up  2 1560.5 1595.2   0.1428                                  
## fit.simple     3 1663.1 1694.4 104.7821     104.64       1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From the ANOVA table we can see the difference in AIC = 102.64 is clearly in favor for the updated simple.up model (Note: a difference in the AIC of two is considered as informative (Burnham and Anderson 2002)).

Further, the difference in the \(\chi^2\) is 104.64–far beyond the threshold value of 3.84 that is necessary to detect a statistically significant difference on a confidence level of \(\alpha = 0.05\).

Saturated model

A saturated model includes all possible paths. As a result, there are no degrees of freedom left and it thus has perfect fit. They represent a special class of model because they allow for everything to add up, meaning we can completely recover the observed matrix of covariances. Unsaturated models have testable implications, however. Under global estimation, our comparison for calculating the GOF is to a saturated model. The reason is that saturated models permit every covariance to be explained, so the model fit function value \(F_{ML}\) fit function goes to zero. Thus, also \(\chi^2 = 0\). In comparison, an unsaturated model will have a positive \(F_{ML}\) (James B. Grace 2021).

Adding each a path from disk to richness (rich) and evenness (even) turns our simple model into a saturated one.

satur <-
"mass.above ~ nadd + disk + rich + even
rich ~ nadd + disk
even ~ nadd + disk

rich ~~ even"

fit.satur <- sem(satur, data = seabloom)
summary(fit.satur, rsq = TRUE)
## lavaan 0.6-9 ended normally after 51 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
##                                                       
##   Number of observations                           240
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   mass.above ~                                        
##     nadd              0.067    0.008    8.118    0.000
##     disk              0.123    0.131    0.941    0.347
##     rich              0.127    0.036    3.523    0.000
##     even              0.361    0.401    0.900    0.368
##   rich ~                                              
##     nadd             -0.129    0.016   -7.999    0.000
##     disk             -0.052    0.291   -0.179    0.858
##   even ~                                              
##     nadd              0.003    0.001    2.363    0.018
##     disk              0.010    0.026    0.374    0.708
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##  .rich ~~                                             
##    .even             -0.261    0.033   -7.916    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .mass.above        0.987    0.090   10.954    0.000
##    .rich              4.871    0.445   10.954    0.000
##    .even              0.040    0.004   10.954    0.000
## 
## R-Square:
##                    Estimate
##     mass.above        0.218
##     rich              0.211
##     even              0.023

Exercise: What modification indices do you expect from a saturated model?

# modindices(fit.satur, minimum.value = 0)

Importance of included paths

In combination with the AIC, from the saturated model it is possible to evaluate if there are any unsupported links in the model.

library("AICcmodavg")

aictab(list(fit.simple, fit.simple.up, fit.satur),
       c("simple", "simple.up", "saturated"))
## 
## Model selection based on AICc:
## 
##            K    AICc Delta_AICc AICcWt Cum.Wt      LL
## simple.up 10 1561.41       0.00   0.89   0.89 -770.22
## saturated 12 1565.68       4.27   0.11   1.00 -770.15
## simple     9 1663.87     102.46   0.00   1.00 -822.54

??? Not sure how this shall show unsupported links…

Summary output

Now, let’s have a closer look on the model via summary requesting also the standardized estimates:

summary(fit.simple.up, fit.measures = TRUE, rsq = TRUE, standardized = TRUE)
## lavaan 0.6-9 ended normally after 45 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        10
##                                                       
##   Number of observations                           240
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 0.143
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.931
## 
## Model Test Baseline Model:
## 
##   Test statistic                               226.074
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000
##   Tucker-Lewis Index (TLI)                       1.039
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -770.223
##   Loglikelihood unrestricted model (H1)       -770.151
##                                                       
##   Akaike (AIC)                                1560.445
##   Bayesian (BIC)                              1595.252
##   Sample-size adjusted Bayesian (BIC)         1563.554
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.037
##   P-value RMSEA <= 0.05                          0.961
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.007
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   mass.above ~                                                          
##     nadd              0.067    0.008    8.118    0.000    0.067    0.529
##     disk              0.123    0.131    0.941    0.346    0.123    0.054
##     rich              0.127    0.036    3.523    0.000    0.127    0.281
##     even              0.361    0.401    0.900    0.368    0.361    0.065
##   rich ~                                                                
##     nadd             -0.129    0.016   -7.998    0.000   -0.129   -0.459
##   even ~                                                                
##     nadd              0.003    0.001    2.362    0.018    0.003    0.151
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .rich ~~                                                               
##    .even             -0.261    0.033   -7.916    0.000   -0.261   -0.594
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .mass.above        0.987    0.090   10.954    0.000    0.987    0.782
##    .rich              4.872    0.445   10.954    0.000    4.872    0.790
##    .even              0.040    0.004   10.954    0.000    0.040    0.977
## 
## R-Square:
##                    Estimate
##     mass.above        0.218
##     rich              0.210
##     even              0.023

Exercise: how would you interpret the meaning of the \(p\)-values of the single paths?

The meaning of the summary output is (James B. Grace 2021):

  • Degrees of freedom represents the number of paths omitted from the model, which provide a capacity to test the architecture of the model
  • p-value: probability of the data given our model
  • Std.err, the standard error
  • Z-value, the analogue to \(t\)-values derived from maximum likelihood estimation (ML)
  • P(>|z|), the \(p\)-values, the probability of obtaining a \(z\) of the given value by chance
  • Regressions = the path coefficients
    • Estimates, the raw unstandardized coefficients
  • Variances = explained variance for endogenous variables = estimates of the error variances
    • Estimates, estimates of the error variances
  • R-square is the variance explained for endogenous variables, the \(1 - error\) variance in standardized terms. Paths from error variables represent influences from un-modeled factors.

The inspect function retrieves information from a lavaan object (for more options see the help):

inspect(fit.simple.up, what = "r2")
## mass.above       rich       even 
##      0.218      0.210      0.023

Visualize the results

The package lavaanPlot allows to simply and straight-forwardly visualize diagrams from lavaan objects. Since it was removed from CRAN begin of 2021, it needs to be downloaded and installed from the archive. The code below first tests, if lavaanPlot is missing in the local library and if so, installs it from the latest archived repository on CRAN.

Another library that handles lavaan objects would be semPlot.

suppressMessages(lavaanPlot_installed <- require(lavaanPlot))
if (!lavaanPlot_installed) {
  install.packages("https://cran.r-project.org/src/contrib/Archive/lavaanPlot/lavaanPlot_0.6.0.tar.gz",
                   repos = NULL, type = "source")
    }

Then, we can plot the results with significance levels displayed as asterisks.

library("lavaanPlot")

lavaanPlot(model = fit.simple.up,
           node_options = list(shape = "box", color = "gray",
                               fontname = "Helvetica"),
           edge_options = list(color = "black"),
           coefs = TRUE, covs = TRUE, stars = c("covs", "regress"))

What is missing in this plot are the \(R^2\) values for the endogenous variables. They can either be represented as the \(R^2\), the variance explained for endogenous variables or as the quantity of error variances (\(\zeta\)) either in raw or standardized units. Another option, if we wish to treat error variables like true causal influences, then we might use path coefficients for their effects. These are the square roots of the error variances (e.g., \(\sqrt{0.84} = 0.92\)) or alternatively, \(\sqrt{1 - R^2}\) (James B. Grace 2021).

Results to report

To present the results of a SEM, following should be stated in the report (James B. Grace 2021):

  • Absolute GOF statistics for final model: \(\chi²\) with the affiliated \(p\)-value, CFI and degrees of freedom (df)
  • Table of raw coefficients and statistics
  • Table of total and indirect effects of interest
  • Computed queries as table or graph

LM vs SEM

The common conclusions for both LM and SEM is that in our example, evenness is of little importance for explaining the total biomass production.

Comparing SEM with LM:

  • SEM represents a network type approach with direct and indirect effects
  • SEM allows for sequential learning
  • LM more likely to detect interactions between variables

References

Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. Second. New York, USA: Springer-Verlag.
Grace, James B. 2006. Structural Equation Modeling and Natural Systems. Cambridge University Press. https://doi.org/10.1017/CBO9780511617799.
Grace, James B. 2021. Quantitative Analysis Using Structural Equation Modeling.” 2021. https://www.usgs.gov/centers/wetland-and-aquatic-research-center/science/quantitative-analysis-using-structural-equation?qt-science_center_objects=0#qt-science_center_objects.
Grace, James B., Donald R. Schoolmaster Jr., Glenn R. Guntenspergen, Amanda M. Little, Brian R. Mitchell, Kathryn M. Miller, and E. William Schweiger. 2012. “Guidelines for a Graph-Theoretic Implementation of Structural Equation Modeling.” Ecosphere 3 (8): art73. https://doi.org/https://doi.org/10.1890/ES12-00048.1.
Hu, Li‐tze, and Peter M. Bentler. 1999. “Cutoff Criteria for Fit Indexes in Covariance Structure Analysis: Conventional Criteria Versus New Alternatives.” Structural Equation Modeling: A Multidisciplinary Journal 6 (1): 1–55. https://doi.org/10.1080/10705519909540118.
Mecklin, Christopher J., and Daniel J. Mundfrom. 2005. “A Monte Carlo Comparison of the Type i and Type II Error Rates of Tests of Multivariate Normality.” Journal of Statistical Computation and Simulation 75 (2): 93–107. https://doi.org/10.1080/0094965042000193233.
Rosseel, Yves. 2012. “Lavaan: An r Package for Structural Equation Modeling.” Journal of Statistical Software 48 (2): 1–36. https://doi.org/10.18637/jss.v048.i02.
———. 2021. lavaan. latent variable analysis.” 2021. https://lavaan.ugent.be/.